/*==============================================================================
案例 5：薪酬影响因素分析（回归算法）
================================================================================

业务场景：
人力资源部门需要了解哪些因素影响员工薪酬，以便制定公平合理的薪酬政策。
通过回归分析识别薪酬的关键驱动因素，并预测不同条件下的合理薪酬水平。

学习目标：
1. 掌握多元线性回归分析
2. 学习 Lasso 回归进行特征选择
3. 理解回归系数的业务含义
4. 掌握回归诊断和模型评估

数据来源：nlsw88.dta（Stata 内置数据集）

作者：Stata ML Course
日期：2025-11-03
==============================================================================*/

clear all
set more off
capture log close

* 设置工作目录
cd "`c(pwd)'"

* 创建输出目录
capture mkdir "output/cases"
capture mkdir "output/cases/figures"
capture mkdir "data/cases"

* 开始日志记录
log using "output/cases/case05_salary_analysis.log", replace text

display "=========================================="
display "案例 5：薪酬影响因素分析"
display "=========================================="
display ""

/*------------------------------------------------------------------------------
第一部分：数据准备
------------------------------------------------------------------------------*/

display "第一部分：数据准备"
display "------------------"

* 加载数据
sysuse nlsw88, clear

* 删除缺失值
drop if missing(wage) | missing(hours) | missing(ttl_exp) | missing(tenure)

display "样本量: " _N
display ""

* 创建新特征
gen monthly_wage = wage * hours * 4  // 月工资估算
gen exp_squared = ttl_exp^2  // 经验平方项（捕捉非线性关系）
gen tenure_squared = tenure^2  // 任期平方项
gen exp_tenure_ratio = ttl_exp / (tenure + 1)  // 经验-任期比率

* 创建分类变量
gen high_education = (grade >= 13)  // 大学及以上
gen experienced = (ttl_exp >= 10)  // 经验丰富
gen long_tenure = (tenure >= 5)  // 长期员工

label variable monthly_wage "月工资（美元）"
label variable exp_squared "工作经验平方"
label variable tenure_squared "任期平方"
label variable exp_tenure_ratio "经验任期比"
label variable high_education "高学历（大学+）"
label variable experienced "经验丰富（10年+）"
label variable long_tenure "长期员工（5年+）"

/*------------------------------------------------------------------------------
第二部分：探索性数据分析
------------------------------------------------------------------------------*/

display ""
display "第二部分：探索性数据分析"
display "------------------------"

* 1. 薪酬分布
summarize wage, detail

histogram wage, ///
    width(1) ///
    frequency ///
    addlabels ///
    title("小时工资分布") ///
    xtitle("小时工资（美元）") ///
    ytitle("频数") ///
    scheme(s2color)
graph export "output/cases/figures/case05_01_wage_distribution.png", replace

* 2. 薪酬与关键变量的关系
twoway (scatter wage ttl_exp, mcolor(blue%30)) ///
       (lfit wage ttl_exp, lcolor(red) lwidth(thick)), ///
    title("工资与工作经验的关系") ///
    xtitle("总工作经验（年）") ///
    ytitle("小时工资（美元）") ///
    legend(order(1 "实际数据" 2 "线性拟合")) ///
    scheme(s2color)
graph export "output/cases/figures/case05_02_wage_vs_experience.png", replace

twoway (scatter wage grade, mcolor(blue%30)) ///
       (lfit wage grade, lcolor(red) lwidth(thick)), ///
    title("工资与教育水平的关系") ///
    xtitle("受教育年限") ///
    ytitle("小时工资（美元）") ///
    legend(order(1 "实际数据" 2 "线性拟合")) ///
    scheme(s2color)
graph export "output/cases/figures/case05_03_wage_vs_education.png", replace

* 3. 不同群体的薪酬对比
graph box wage, over(union) over(married) ///
    title("不同群体的工资分布") ///
    subtitle("按工会状态和婚姻状态分组") ///
    ytitle("小时工资（美元）") ///
    scheme(s2color)
graph export "output/cases/figures/case05_04_wage_by_groups.png", replace

* 4. 相关性分析
correlate wage ttl_exp tenure grade hours age

/*------------------------------------------------------------------------------
第三部分：数据分割
------------------------------------------------------------------------------*/

display ""
display "第三部分：数据分割"
display "------------------"

set seed 20251103
gen random = runiform()
gen train = (random <= 0.8)

count if train
display "训练集样本量: " r(N)
count if !train
display "测试集样本量: " r(N)

/*------------------------------------------------------------------------------
第四部分：模型 1 - 基础 OLS 回归
------------------------------------------------------------------------------*/

display ""
display "第四部分：基础 OLS 回归"
display "----------------------"

* 简单模型（只包含核心变量）
regress wage ttl_exp grade hours if train

* 保存模型
estimates store ols_simple

* 显示结果
display ""
display "简单 OLS 模型结果："
display "-------------------"
estimates replay ols_simple

* 解释关键系数
display ""
display "关键发现："
display "- 每增加 1 年工作经验，工资增加 " %5.3f _b[ttl_exp] " 美元/小时"
display "- 每增加 1 年教育，工资增加 " %5.3f _b[grade] " 美元/小时"
display "- R² = " %5.3f e(r2)

/*------------------------------------------------------------------------------
第五部分：模型 2 - 完整 OLS 回归
------------------------------------------------------------------------------*/

display ""
display "第五部分：完整 OLS 回归"
display "----------------------"

* 包含更多变量和交互项
regress wage ttl_exp exp_squared grade hours tenure married union south ///
    c.ttl_exp#c.grade if train

estimates store ols_full

* 显示结果
display ""
display "完整 OLS 模型结果："
display "-------------------"
estimates replay ols_full

* 模型诊断
predict wage_pred_ols if train
predict resid_ols if train, residuals

* 残差图
twoway (scatter resid_ols wage_pred_ols, mcolor(blue%30)), ///
    title("OLS 回归残差图") ///
    xtitle("预测工资（美元/小时）") ///
    ytitle("残差") ///
    yline(0, lcolor(red) lpattern(dash)) ///
    scheme(s2color)
graph export "output/cases/figures/case05_05_ols_residuals.png", replace

/*------------------------------------------------------------------------------
第六部分：模型 3 - Lasso 回归（特征选择）
------------------------------------------------------------------------------*/

display ""
display "第六部分：Lasso 回归"
display "-------------------"

* Lasso 回归（交叉验证选择 lambda）
lasso linear wage ttl_exp exp_squared grade hours tenure tenure_squared ///
    age married union south industry c.ttl_exp#c.grade ///
    c.ttl_exp#i.married c.grade#i.union ///
    if train, selection(cv) rseed(20251103)

* 显示选择的变量
lassocoef, display(coef, standardized)

* 预测
predict wage_pred_lasso if !train

/*------------------------------------------------------------------------------
第七部分：模型 4 - Elastic Net 回归
------------------------------------------------------------------------------*/

display ""
display "第七部分：Elastic Net 回归"
display "-------------------------"

* Elastic Net（alpha=0.5，平衡 L1 和 L2 正则化）
elasticnet linear wage ttl_exp exp_squared grade hours tenure tenure_squared ///
    age married union south industry c.ttl_exp#c.grade ///
    if train, alpha(0.5) selection(cv) rseed(20251103)

* 预测
predict wage_pred_enet if !train

/*------------------------------------------------------------------------------
第八部分：模型对比与评估
------------------------------------------------------------------------------*/

display ""
display "第八部分：模型对比与评估"
display "------------------------"

* 在测试集上预测（OLS）
quietly regress wage ttl_exp exp_squared grade hours tenure married union south ///
    c.ttl_exp#c.grade if train
predict wage_pred_ols_test if !train

* 计算评估指标
foreach model in ols_test lasso enet {
    quietly gen sq_error_`model' = (wage - wage_pred_`model')^2 if !train
    quietly gen abs_error_`model' = abs(wage - wage_pred_`model') if !train
    
    quietly summarize sq_error_`model'
    local mse_`model' = r(mean)
    local rmse_`model' = sqrt(r(mean))
    
    quietly summarize abs_error_`model'
    local mae_`model' = r(mean)
    
    quietly correlate wage wage_pred_`model' if !train
    local r_`model' = r(rho)
    local r2_`model' = r(rho)^2
}

* 显示对比结果
display ""
display "模型性能对比（测试集）"
display "======================"
display ""
display "模型          RMSE      MAE       R²"
display "--------------------------------------"
display "OLS         " %7.3f `rmse_ols_test' "   " %7.3f `mae_ols_test' "   " %6.3f `r2_ols_test'
display "Lasso       " %7.3f `rmse_lasso' "   " %7.3f `mae_lasso' "   " %6.3f `r2_lasso'
display "Elastic Net " %7.3f `rmse_enet' "   " %7.3f `mae_enet' "   " %6.3f `r2_enet'
display ""

* 预测 vs 实际对比图
twoway (scatter wage wage_pred_ols_test if !train, mcolor(blue%30)) ///
       (function y=x, range(0 30) lcolor(red) lpattern(dash)), ///
    title("OLS 模型：预测 vs 实际工资") ///
    subtitle("测试集") ///
    xtitle("预测工资（美元/小时）") ///
    ytitle("实际工资（美元/小时）") ///
    legend(order(1 "数据点" 2 "完美预测线")) ///
    scheme(s2color)
graph export "output/cases/figures/case05_06_prediction_vs_actual.png", replace

/*------------------------------------------------------------------------------
第九部分：业务洞察与建议
------------------------------------------------------------------------------*/

display ""
display "第九部分：业务洞察与建议"
display "------------------------"

* 1. 关键薪酬驱动因素
display ""
display "1. 关键薪酬驱动因素（基于 OLS 完整模型）"
display "----------------------------------------"

quietly regress wage ttl_exp exp_squared grade hours tenure married union south ///
    c.ttl_exp#c.grade if train

display ""
display "主要发现："
display "- 教育水平：每增加 1 年教育，工资增加 " %5.2f _b[grade] " 美元/小时"
display "- 工作经验：每增加 1 年经验，工资增加 " %5.2f _b[ttl_exp] " 美元/小时"
display "- 工会成员：工会成员工资高出 " %5.2f _b[union] " 美元/小时"
display "- 婚姻状态：已婚员工工资高出 " %5.2f _b[married] " 美元/小时"

* 2. 薪酬预测示例
display ""
display "2. 薪酬预测示例"
display "----------------"

* 场景 1：新员工（本科，无经验）
local pred1 = _b[_cons] + _b[grade]*16 + _b[ttl_exp]*0 + _b[hours]*40
display "场景 1 - 本科新员工（16年教育，0年经验，40小时/周）"
display "   预测工资: " %6.2f `pred1' " 美元/小时"
display "   月工资估算: " %8.2f (`pred1'*40*4) " 美元"

* 场景 2：经验丰富员工（本科，10年经验）
local pred2 = _b[_cons] + _b[grade]*16 + _b[ttl_exp]*10 + _b[exp_squared]*100 + _b[hours]*40
display ""
display "场景 2 - 经验丰富员工（16年教育，10年经验，40小时/周）"
display "   预测工资: " %6.2f `pred2' " 美元/小时"
display "   月工资估算: " %8.2f (`pred2'*40*4) " 美元"

* 3. 薪酬公平性分析
display ""
display "3. 薪酬公平性分析"
display "----------------"

* 计算残差（实际工资 - 预测工资）
predict wage_expected if train, xb
gen wage_gap = wage - wage_expected if train

* 识别薪酬过高/过低的员工
gen overpaid = (wage_gap > 2) if train
gen underpaid = (wage_gap < -2) if train

count if overpaid
display "薪酬过高员工数量: " r(N)
count if underpaid
display "薪酬过低员工数量: " r(N)

/*------------------------------------------------------------------------------
第十部分：保存结果
------------------------------------------------------------------------------*/

display ""
display "第十部分：保存结果"
display "------------------"

* 保存分析结果
preserve
keep if !train
keep wage wage_pred_ols_test wage_pred_lasso wage_pred_enet
export delimited using "output/cases/case05_predictions.csv", replace
display "预测结果已保存: output/cases/case05_predictions.csv"
restore

/*------------------------------------------------------------------------------
总结
------------------------------------------------------------------------------*/

display ""
display "=========================================="
display "案例 5 完成！"
display "=========================================="
display ""
display "主要成果："
display "1. 完成薪酬影响因素的多元回归分析"
display "2. 对比 OLS、Lasso、Elastic Net 三种模型"
display "3. 识别关键薪酬驱动因素"
display "4. 提供薪酬预测和公平性分析"
display "5. 生成 6 个可视化图表"
display ""
display "关键发现："
display "- 教育和经验是薪酬的主要驱动因素"
display "- 工会成员和已婚员工薪酬更高"
display "- 模型 R² 约为 " %5.3f `r2_ols_test' "，具有较好的预测能力"
display ""
display "输出文件："
display "- 图表: output/cases/figures/case05_*.png (6个)"
display "- 预测: output/cases/case05_predictions.csv"
display "- 日志: output/cases/case05_salary_analysis.log"
display ""

log close

